\documentclass[11pt]{beamer}
\usepackage[UTF8,scheme=plain]{ctex}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usetheme{Boadilla}

\usepackage{framed} % 可以用 \begin{shaded}，即背景色块
\definecolor{shadecolor}{rgb}{0.9,0.9,0.9}

\lstset{
	basicstyle=\footnotesize,
    flexiblecolumns,
    showstringspaces = false,
	numbers=left,
    numberstyle=\tiny,
    numbersep=4pt,
	breaklines,%自动换行
	frame=shadowbox,%加边框
	backgroundcolor = \color{red!30!green!30!blue!30},
	%背景色：蓝灰色
	language=C++, % 语言为 c++
	commentstyle=\color{red!80!green!80!blue!80},% 注释用蓝灰色字体标出
	keywordstyle=\color{blue!70}, % 高亮部分用蓝色字体
	tabsize = 4,
    escapeinside=@@
}

\begin{document}
	\author{ 路毅 \hspace{0.3cm} 曲阜师范大学 }
	\date{\number\year 年 \number\month 月 \number\day 日}
	\title{计算物理第4部分：微分方程}
	%\subtitle{}
	%\logo{}
	%\institute{}
	%\date{}
	%\subject{}
	%\setbeamercovered{transparent}
	%\setbeamertemplate{navigation symbols}{}
	\begin{frame}[plain]
	\maketitle
\end{frame}

\kaishu

\begin{frame}[fragile]{微分方程}
常微分方程，即一元函数的微分方程，比如自由落体方程
\begin{equation}
\frac{d^2}{dt^2} x = - g,
\end{equation}
偏微分方程，就是多元函数的微分方程，比如静电势方程
\begin{equation}
\nabla^2 \psi(x,y,z) = - \frac{\rho(x,y,z)}{\epsilon_0},
\end{equation}
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十章：常微分方程 }
            \begin{itemize}
            \item {\color{blue} 初值问题}
            \item 边值问题
            \end{itemize}
		\item 第十一章：偏微分方程
	\end{itemize}
\end{frame}

\begin{frame}{初值问题}
如果函数$y$在定义域一端$x=a$处的信息已知，需要推测其他地方的函数值，那么这种问题就叫做初值问题。
一个典型的例子是，将单摆拉开一定角度，然后松手，推测一天以后单摆在哪里，速度是多少。
将松手的瞬间定位$t=0$，则
\begin{eqnarray}
\theta(t=0) &=& \theta_0, \\
\dot{\theta}(t=0) &=& 0
\end{eqnarray}
定义了初始时刻的摆角和角速度。
然后根据牛顿第二定律做动力学的分析，推测$\theta(t)$与$\dot{\theta}(t)$的轨迹，这就是初值问题。
\end{frame}

\begin{frame}{边值问题}
相对的，边值问题是这样的一类问题，我们知道定义域一端$x=a$处$y$的一部分信息，以及另一端$x=b$的一部分信息，然后根据动力学方程推断出整个定义域内各处的$y$值。
可怜的中学生们也可能遇到这样的问题，比如，$t=0$时刻，单摆的摆角为
\begin{equation}
\theta(t=0) = \theta_1,
\end{equation}
但是不知道它的角速度是多少。
而在$t=100s$，单摆的摆角为
\begin{equation}
\theta(t=10) = \theta_2,
\end{equation}
也不知道角速度是多少，求单摆的最大摆角$\theta_0$。
这就是一个时间轴上的边值问题。
\end{frame}

\begin{frame}{微分方程降阶}
如果方程中最高阶的导数是$n$阶导数$y^{(n)}$，我们就把这个方程叫做$n$阶微分方程，
\begin{equation}
\frac{d^n}{dx^n} y = f(y, \frac{dy}{dx}, \cdots, \frac{d^{n-1}}{dx^{n-1}}y),
\end{equation}
一般需要$n-1$个关于$y,y',y'',\cdots$的初始条件。
为了方便处理，我们一般把高阶的微分方程改写为一阶微分方程组，定义$n$个函数$y^{0}, \cdots, y^{(n-1)}$，分别表示函数$y$的$0,\cdots,n-1$阶导函数。
\end{frame}

\begin{frame}{微分方程降阶}
则原来的微分方程可以写作，
\begin{eqnarray}
\frac{d}{dx} y^{(0)} &=& y^{(1)}, \\
\frac{d}{dx} y^{(1)} &=& y^{(2)}, \\
& \vdots &  \\
\frac{d}{dx} y^{(n-1)} &=& f(y^0, \cdots, y^{(n-1)}), \\
\end{eqnarray}
这$n$个1阶微分方程，加上$n-1$个关于$y^{(0)},y^{(1)},\cdots,y^{(n-1)}$的初始条件，唯一地确定了解。
\end{frame}

\begin{frame}{微分方程降阶}
为了简单形象，我们举一个例子，一维谐振子，即光滑平面上固定在横向理想弹簧上的滑块，本来在平衡位置，现在拉开一段距离$a$（弹性形变），然后松开，谐振子的运动方程为
\begin{equation}
m \frac{d^2}{dt^2} x(t) = - kx(t),
\end{equation}
为了简便，我们设$m=k=a=1$（也可以通过定义新的变量将方程变为无量纲方程，效果是一样，简便）。
那么方程变为
\begin{equation}
\frac{d^2}{dt^2} x = -x,
\end{equation}
初始条件为$x(0)=1, x'(0)=0$，为了降阶，我们定义$v(t)=\frac{dx}{dt}$，方程变为
\begin{eqnarray}
\frac{dx}{dt} &=& v, ~~ x(0)=1, \\
\frac{dv}{dt} &=& - x, ~~ v(0)=0.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{初值问题}
我们的目标是求得$y(x)$在定义域中每一点的值。
但有时没有解析解，于是我们退一步求其次，在定义域中取一系列离散的点，近似地求得这些点上的函数值，也可以绘制出函数的样子，得到对函数的了解。
\vspace{0.5cm}
所以我们将定义域$\left[x_0=a, x_n=b\right]$等分为$n$份，步长
\begin{equation}
h=(b-a)/n,
\end{equation}
得到离散点
\begin{equation}
x_0, x_1, x_2, \cdots, x_{n-1}, x_n,
\end{equation}
然后求取$f(x_i)$的近似值$y_i = f(x_i)$，就算是完成了数值任务——计算机高兴地说“我解出来了！”。
\end{frame}

\begin{frame}[fragile]{欧拉法}
最简单直接的是欧拉法，做$y(x)$的一阶泰勒展开，
\begin{eqnarray}
y(x_{i+1}) = y(x_i) + h \frac{dy}{dx}|_{x_i} + O(h^2) = y(x_i)+ hf(x_i, y(x_i)) + O(h^2),
\end{eqnarray}
所以可以构造迭代公式
\begin{eqnarray}
y_{i+1} = y_i + hf(x_i, y_i)
\end{eqnarray}
完成 $y_i$ 到 $y_{i+1}$ 的迭代。
注意，上面我们用$y(x_i)$代表函数值的理论值，用$y_i$代表$y(x_i)$的数值近似。
在每一步迭代，截断误差约为 $O(h^2)$，所以迭代求得 $y_n$ 时，误差为
\begin{equation}
y_n - y(x_n) \propto nh^2 \propto h,
\end{equation}
所以欧拉法的精度为$O(h)$，即一阶精度。
\end{frame}

\begin{frame}[fragile]{欧拉法}
我们可以想到一个很简单的例子，即雨滴下落。
雨滴从云中掉落，重力使它持续加速，如果没有空气阻力，就会像水做的子弹一样，把人类砸伤。

还好有空气阻力，空气阻力与雨滴速度的平方成正比，使得雨滴速度增加到某个值，就不再增加了。
所以雨滴落地时速度一般是几米每秒，人类可以在雨中愉快地歌唱。

\begin{figure}
\includegraphics[height=0.5\textheight]{singing-in-the-rain}
\end{figure}
\end{frame}

\begin{frame}[fragile]{欧拉法}
雨滴速度的初始值是
\begin{equation}
v(0) = 0,
\end{equation}
其微分方程为
\begin{equation}
\frac{dv}{dt} = g - \alpha v^2,
\end{equation}
其中$\alpha$是参数，与雨滴质量等因素有关，我们简单将它设为1，也许不准确，但不妨碍我们展示这个小故事。
\end{frame}

\begin{frame}[fragile]{欧拉法}
我们用欧拉法解微分方程，确定$v(t)$，代码如下
\lstset{basicstyle=\footnotesize}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

/* Euler(...)           solve 1-order ODE for initial problem
 * double a             [a,b]
 * double b
 * double *func(double x, double y)     f(x,y)=dy/dx
 * int n                [a,b] is divided into n pieces
 * double y0            initial value y(a)
 * double *y            y[n+1] stores y0, y1,..., yn
 */
void Euler(double a, double b, double (*func)(double x, double y), int n, double y0, double *y){
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{欧拉法} 
\lstset{firstnumber = last}   
\begin{lstlisting}   
        int i;
        double h = (b-a)/n; //步长
        y[0] = y0;
        for(i=0;i<n;i++){
                y[i+1] = y[i] + h*func(a + i*h, y[i]);
        }
        ofstream fp("Euler.txt"); //将结果输出到Euler.txt，用作画图
        for(i=0;i<=n;i++){
                fp<<a+i*h<<"\t"<<y[i]<<endl;
        }
        fp.close();
}

double f(double x, double y){ // dv/dt = 10 - v^2
        return 10 - y*y;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{欧拉法} 
\lstset{firstnumber = last}   
\begin{lstlisting}  
int main(){
        
        int i,j,k;
        
        int n=10000;
        double *y = new double [n];
        Euler(0, 2, f, n, 0, y);
        
        delete [] y;
        return 0;
} 
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{欧拉法} 
跑完程序以后，我们得到程序文件夹下的Euler.txt文件，就可以用gnuplot画图，gnuplot定位到当前文件夹后，打入
\begin{lstlisting}
plot "Euler.txt" w l
\end{lstlisting}
即用折线连接Euler.txt文件中给出的离散点，下图中我们可以看到，雨滴速度变大，然后趋于饱和的过程，最终速度为3m/s左右，所以是一场小雨。
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Euler}
%\caption{雨滴速度(m/s) v.s. 时间(s)。}
%\label{fig:Euler-raindrop}
\end{figure}
\end{frame}

\begin{frame}[fragile]{欧拉法} 
我们可以换一个角度看这个问题，对上面描述的这种初值问题，$y(x_{i+1})$的积分形式为
\begin{equation}
y(x_{i+1}) = y(x_i) + \int_{x_i}^{x_{i+1}} f(x,y) dx,
\end{equation}
上式右边的积分，如果做如下矩形近似
\begin{eqnarray}
\int_{x_i}^{x_{i+1}} f(x,y) dx &=& \int_{x_i}^{x_{i+1}} (f(x_i,y(x_i)) + O(x) ) dx   \nonumber \\
&\approx& h f(x_i,y(x_i)) + O(h^2),
\end{eqnarray}
则得到欧拉法。
\end{frame}

\begin{frame}[fragile]{改进的欧拉法} 
那么如果我们希望得到更好的近似，可以用梯形法构造上述积分，
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = \frac{h}{2}[ f(x_i,y(x_i)) + f(x_{i+1}, y(x_{i+1}) ) ] + O(h^3),
\end{equation}
但是我们不知道$y(x_{i+1})$，所以做近似
\begin{equation}
y(x_{i+1}) = y(x_i) + h f(x_i,y_i) + O(h^2),
\end{equation}
得到
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = \frac{h}{2}[ f(x_i,y(x_i)) + f(x_{i+1}, y(x_i) + h f(x_i,y_i)) ]  + O(h^3),
\end{equation}
用相应的数值近似代替上式中的解析表达式，得到迭代公式
\begin{equation}
y_{i+1} = y_i + \frac{h}{2}[ f(x_i,y_i) + f(x_{i+1}, y_i + h f(x_i,y_i)) ],
\end{equation}
这就是改进的欧拉法，每步具有$O(h^3)$精度，推到$y_n$时，具有$O(h^2)$精度。
\end{frame}

\begin{frame}[fragile]{龙格库塔方法} 
延续上面的思路，如果用中值法处理积分，则得到二阶龙格库塔方法，
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = h f(x_i + h/2, y(x_i+h/2)) + O(h^3),
\end{equation}
未知的$f(x_i + h/2, y(x_i+h/2))$可以用下面的近似公式替代，
\begin{eqnarray}
f(x_i + \frac{h}{2}, y(x_i+h/2)) &=& f(x_i+\frac{h}{2}, y(x_i) + \frac{h}{2}f(x_i,y(x_i)) + O(h^2)) \nonumber\\
&=& f(x_i+\frac{h}{2}, y(x_i) + \frac{h}{2}f(x_i,y(x_i))) + O(h^2),
\end{eqnarray}
得到公式
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = h f(x_i + \frac{h}{2}, y(x_i) + \frac{h}{2}f(x_i, y(x_i)) ) + O(h^3),
\end{equation}
\end{frame}

\begin{frame}[fragile]{龙格库塔方法} 
相应的，数值迭代公式为
\begin{equation}
y_{i+1} = y_i + h f(x_i + \frac{h}{2}, y_i + \frac{h}{2} f(x_i,y_i)),
\end{equation}
或者写得稍微清晰一些
\begin{equation}
\left\{
\begin{aligned}
&K_1 = f(x_i, y_i)\\
&y_{i+1} = y_i + h f(x_i + h/2, y_i + (h/2)K_1)
\end{aligned}
\right.
\end{equation}
这个公式每步具有$O(h^3)$，推到$y_n$时累积有$O(h^2)$精度，所以叫做二阶龙格库塔公式。
\end{frame}

\begin{frame}[fragile]{4阶龙格库塔} 
若需要四阶龙格库塔公式，即需要每步精度为$O(h^5)$，我们想起了辛普生方法，即
\begin{eqnarray}
\int^{x_i +h}_{x_i} f(x,y) dx &=& \frac{h}{6}[ f(x_i,y(x_i)) + 4f(x_i + \frac{h}{2},y(x_i+\frac{h}{2}))  \nonumber\\
&& + f(x_i + h,y(x_i+h)) ] + O(h^5),
\end{eqnarray}
找到$f(x_i + \frac{h}{2},y(x_i+\frac{h}{2}))$，$f(x_i + h,y(x_i+h))$的4阶精度替代方程，即可构造出4阶龙格库塔公式，下面的数值迭代式是前人的构造，
\begin{eqnarray}
\left\{
\begin{aligned}
&K_1 = f(x_i, y_i)\\
&K_2 = f(x_i + h/2, y_i + (h/2)K_1)\\
&K_3 = f(x_i + h/2, y_i + (h/2)K_2)\\
&K_4 = f(x_i + h, y_i + hK_3)\\
&y_{i+1} = y_i + \frac{h}{6}(K_1 + 2K_2 + 2K_3 + K_4)
\end{aligned}
\right.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{4阶龙格库塔}
如果是2阶以上的常微分方程，则降阶以后会得到一系列1阶常微分方程。
不妨设$\vec{y} = \{ y^{(0)}, y^{(1)}, \cdots, y^{(n-1)} \}^\top $，假设微分方程为
\begin{eqnarray}
\frac{d}{dx} y^{(0)} &=& y^{(1)} = f_0(\vec{y}), \\
\frac{d}{dx} y^{(1)} &=& y^{(2)} = f_1(\vec{y}), \\
\frac{d}{dx} y^{(2)} &=& y^{(3)} = f_2(\vec{y}), \\
                    &\vdots& \\
\frac{d}{dx} y^{(n-1)} &=& y^{(n)} = f_{n-1}(\vec{y}).
\end{eqnarray}
即
\begin{equation}
\frac{d}{dx} \vec{y} = \vec{f}(\vec{y}).
\end{equation}
\end{frame}

\begin{frame}[fragile]{4阶龙格库塔}
其中
\begin{equation}
\vec{f}(\vec{y}) = \left\{ f_0(\vec{y}), f_1(\vec{y}), \cdots, f_{n-1}(\vec{y}) \right\}^\top
\end{equation}
高阶常微分方程的4阶龙格库塔方法可以写作
\begin{eqnarray}
\left\{
\begin{aligned}
&\vec{K}_1 = \vec{f}(x_i, \vec{y}_i)\\
&\vec{K}_2 = \vec{f}(x_i + h/2, \vec{y}_i + (h/2)\vec{K}_1)\\
&\vec{K}_3 = f(x_i + h/2, \vec{y}_i + (h/2)\vec{K}_2)\\
&\vec{K}_4 = f(x_i + h, \vec{y}_i + h\vec{K}_3)\\
&\vec{y}_{i+1} = \vec{y}_i + \frac{h}{6}(\vec{K}_1 + 2\vec{K}_2 + 2\vec{K}_3 + \vec{K}_4)
\end{aligned}
\right.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{4阶龙格库塔方法} 
这个公式在每个步长的误差为$O(h^5)$，累积误差为$O(h^4)$。
上面这个式子中，$K1$是$x_i$处的斜率，而$K2,K3$是$x_i+h/2$处斜率的两个近似，$K4$是$x_i+h$处的斜率。
所以四阶龙格库塔方法利用了3处的斜率信息，然后进行加权，所以能得到对$y_{i+1}$的更精确的推断。
四阶龙格库塔的应用相当普遍，下面的代码实现一个步长的计算
\begin{lstlisting}
/* rk4:                 4阶龙格库塔一个步长的计算
 * double x:            xi
 * double h:            步长
 * double *y:           xi处已知的 y, y', y'', ..., y^{(n-1)}
 * double *ynext        计算的结果 xi+h 处的 y, y', y'', ..., y^{(n-1)}
 * void (*deriv)(int n, double x, double *y, double *dydx)
 *                      xi 处的导函数，根据微分方程分别计算 y, y', ..., y^{(n-1)} 在 x 处的导函数，并
储存在 double *dydx 中
 */
void rk4(double x, double h, int n, double *y, double *ynext, void (*deriv)(int n, double x, double *y, double *dydx)){
        int i,j;
        double hh = h*0.5;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{龙格库塔方法} 
\lstset{firstnumber = last}
\begin{lstlisting}
        double h6 = h/6;
        double *K, *ytemp;
        K = new double [n];
        deriv(n, x, y, K);// K1
        for(i=0;i<n;i++) ynext[i] = K[i];

        ytemp = new double [n];
        for(i=0;i<n;i++) ytemp[i] = y[i] + hh*K[i];// yi + h/2*K1
        deriv(n, x+hh, ytemp, K);// K2 
        for(i=0;i<n;i++) ynext[i] += 2 * K[i];

        for(i=0;i<n;i++) ytemp[i] = y[i] + hh*K[i];// yi + h/2*K2
        deriv(n, x+hh, ytemp, K);// K3 
        for(i=0;i<n;i++) ynext[i] += 2 * K[i];

        for(i=0;i<n;i++) ytemp[i] = y[i] + h*K[i];// yi + h*K3
        deriv(n, x+h, ytemp, K);// K4
        for(i=0;i<n;i++) ynext[i] += K[i];
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{龙格库塔方法} 
\lstset{firstnumber = last}
\begin{lstlisting}

        for(i=0;i<n;i++){// update y_{j+1}
                ynext[i] = y[i] + h6*ynext[i];
        }
        delete [] K;
        delete [] ytemp;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{经典谐振子} 
下面演示怎样用4阶龙格库塔方法，求解经典谐振子，前面已经提到，谐振子的运动方程为
\begin{equation}
m \frac{d^2}{dt^2} x(t) = - kx(t),
\end{equation}
为了简便，我们设$m=k=a=1$，则方程变为
\begin{equation}
\frac{d^2}{dt^2} x = -x,
\end{equation}
初始条件为$x(0)=1, x'(0)=0$，则有
\begin{eqnarray}
\frac{dx}{dt} &=& v, \\
\frac{dv}{dt} &=& - x.
\end{eqnarray}
然后利用龙格库塔方法，即可实现求解，代码如下
\end{frame}

\begin{frame}[fragile]{经典谐振子}
\begin{lstlisting}
头文件声明。。。

void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];// dx/dt = v
        dydx[1] = - y[0];// dv/dt = -x
}

int main(){
        int i,j, n=2, N=1000;
        double a=0, b=10;
        double h = (b-a)/N; //步长
        double *y = new double [n]; //每步已知的 x,v
        double *ynext = new double [n]; //每步未知的 x,v
        ofstream fp("classic_HO.txt");//输出文件流对象，将每步求得的结果输出到文件
        y[0]=1; y[1] = 0;//初始值 x(0)=1, v(0) = 0
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{龙格库塔方法} 
\lstset{firstnumber = last}
\begin{lstlisting}
        fp<<a<<"\t"<<y[0]<<"\t"<<y[1]<<"\n"; // 输出t=0时的x,v
        for(i=0;i<N;i++){
                rk4(a+i*h, h, n, y, ynext, derivative);
                for(j=0;j<n;j++) y[j] = ynext[j];
                fp<<a+(i+1)*h<<"\t"<<ynext[0]<<"\t"<<ynext[1]<<endl;
                //把每步的结果，即 a+(i+1)*h 处的 x,v 输出
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{经典谐振子}
得到的数据存放在"classic\_Ho.txt"中，然后可以用gnuplot画图
\begin{lstlisting}
plot "classic_HO.txt" u 1:2 w l lw 3 t "x (m), Arial, 15", "" u 1:3 w l lw 3 t "v (m/s)"
\end{lstlisting}
上面的gnuplot指令表示，利用文件中的第1、2列分别作为横、纵坐标画一条曲线，第1、3列分别作为横、纵坐标画另一条曲线。得到如下图片
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Classic_HO}
\label{fig:ClassicHO}
%\caption{经典谐振子：4阶龙格库塔求解。}
\end{figure}
\end{frame}

\begin{frame}{阻尼振动}
如果谐振子浸泡在液体中，它的运动就会受到液体的阻力，运动速度较小时，阻力与它的速度成正比，方向相反，
\begin{equation}
F_f = - \gamma \frac{dx}{dt} = - \gamma v,
\end{equation}
所以振子的运动方程变为
\begin{equation}
m \frac{d^2}{dt^2} x = - k x - \gamma v,
\end{equation}
因为我们假设$m = k = 1$，所以方程变为
\begin{equation}
\frac{d^2}{dt^2} x = - x - \gamma v,
\end{equation}
方程降阶以后得到
\begin{eqnarray}
\frac{d}{dt} x &=& v, \\
\frac{d}{dt} v &=& -x - \gamma v.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{阻尼振动}
我们只需要稍微修改上面代码的$derivative$函数，即可计算阻尼运动。
\begin{lstlisting}
void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];// dx/dt = v
        dydx[1] = - y[0] - 0.1 * y[1];// dv/dt = -x - gamma v
}
\end{lstlisting}
为了看得更清楚，我们计算$t=[0,100]s$的振动，得到图中所示的结果。
\begin{lstlisting}
set term png
set output "classicHO_damping.png"
set xlabel "t(s)" font "Arial,15"
set title "Damping of classic harmonic oscillator" font "Arial,15"
plot "classicHO_damping.txt" u 1:2 w l lw 3 t "x(m)"#, "" u 1:3 w l lw 3 t "v(m/s)"
\end{lstlisting}
\end{frame}

\begin{frame}{阻尼振动}
可以看到，因为液体阻尼的存在，简谐振动的振幅快速衰减。
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{classicHO_damping}
\label{fig:HOdamping}
\caption{谐振子阻尼运动，$m=k=a=1$，$\gamma=0.1$。}
\end{figure}
\end{frame}

\begin{frame}{共振}
如果除了阻尼以外，谐振子还受一个周期性的驱动力，
\begin{equation}
F = F_0 \sin(\omega t),
\end{equation}
那么谐振子的运动方程就变为
\begin{equation}
m \frac{d^2}{dt^2} x = - k x - \gamma v + F_0 \sin(\omega t),
\end{equation}
如果我们假设$m=k=a=1$，则有
\begin{eqnarray}
\frac{d}{dt} x &=& v, \\
\frac{d}{dt} v &=& -x - \gamma v + F_0 \sin(\omega t).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{共振}
我们只需要稍微修改上面代码的$derivative$函数，即可计算阻尼运动。
下面是$\gamma = 0.1, F_0 = 0.05$的例子，
\begin{lstlisting}
void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];// dx/dt = v
        dydx[1] = - y[0] - 0.1 * y[1] + 0.05 * sin(1.1*x);// dv/dt = -x - gamma v + F_0 sin(omega t)
}
\end{lstlisting}
\end{frame}


\begin{frame}
我们设$\gamma=0.1, F_0 = 0.05$，分别取$\omega=0.8,0.9,1.0,1.1$，得到下图。
\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{classicHO_Resonance}
\label{fig:classicHO_Resonance}
%\caption{谐振子在外部驱动力下的共振现象。}
\end{figure}
\end{frame}

\begin{frame}{Numerov方法}
对于如下形式的常微分方程
\begin{eqnarray}
 \frac{ d^2 y}{dx^2} = - g(x) y(x) + s(x),
 \label{Numerov1}
\end{eqnarray}
其中 $g(x), s(x)$ 的形式是已知的，前苏联天文学家 Numerov 提出了每步有$O(h^6)$精度的算法：
\begin{eqnarray}
&&y(x_{n+1}) [ 1 + \frac{h^2}{12} g(x_{n+1}) ] \nonumber\\
&=& 2 y(x_n) [ 1 - \frac{ 5 h^2}{12} g(x_n) ] - y(x_{n-1}) [1 + \frac{h^2}{12} g(x_{n-1}) ]
\nonumber\\
&& + \frac{ h^2 }{12} [ s(x_{n+1}) + 10 s(x_n) + s(x_{n-1}) ] + O(h^6),
\end{eqnarray}
如果用节点上的计算值 $y_{n+1}, y_n, y_{n-1}$ 代替上式中的 $y(x_{n+1}), y(x_n), y(x_{n-1})$，并忽略上式右边最后一项 $O(h^6)$，则得到一个精度为$O(h^6)$的数值递推式
\end{frame}

\begin{frame}{Numerov方法}
\begin{eqnarray}
y_{n+1} ( 1 + \frac{h^2}{12} g_{n+1}) &=& 2 y_n ( 1 - \frac{ 5 h^2}{12} g_n ) - y_{n-1} (1 + \frac{h^2}{12} g_{n-1} ) \nonumber\\
&& + \frac{ h^2 }{12} ( s_{n+1} + 10 s_n + s_{n-1} ),
\end{eqnarray}
其中 $g_n = g(x_n), ~ s_n = s(x_n)$。
如果区间最左边的两个节点函数值 $y_0, y_1$ 已知，就可以推得 $y_2, y_3, \cdots, y_n$，即得到方程的数值解。
\end{frame}

\begin{frame}{Numerov方法}
下面我们推导 Numerov 方法的公式。
$y(x)$ 的泰勒展开式有
\begin{eqnarray}
y( x_{n+1} ) &=& y(x_n) + h y'(x_n) + \frac{ h^2}{2} y''(x_n) + \frac{ h^3}{6} y'''(x_n) 
\nonumber\\
&& + \frac{ h^4}{24} y''''(x_n) + \frac{ h^5}{120} y'''''(x_n) + O( h^6),
\nonumber\\
y( x_{n-1} ) &=& y(x_n) - h y'(x_n) + \frac{ h^2}{2} y''(x_n) - \frac{ h^3}{6} y'''(x_n) 
\nonumber\\
&& + \frac{ h^4}{24} y''''(x_n) - \frac{ h^5}{120} y'''''(x_n) + O( h^6),
\nonumber\\
\end{eqnarray}
上两式相加，得到
\begin{eqnarray}
y( x_{n+1} ) + y( x_{n-1} ) - 2 y(x_n) = h^2 y'' (x_n) + \frac{ h^4 }{12} y''''(x_n) + O(h^6).
\label{Numerov2}
\end{eqnarray}
\end{frame}

\begin{frame}{Numerov方法}
结合（\ref{Numerov1}），可以推得
\begin{eqnarray}
&& y( x_{n+1} ) + y( x_{n-1} ) - 2 y(x_n) \nonumber\\
&=& h^2 ( - g(x_n) y(x_n) + s(x_n) ) + \frac{ h^4 }{12} y''''(x_n) + O(h^6).
\label{Numerov3}
\end{eqnarray}
$ y''''(x_n) $ 可以用（\ref{Numerov2}）类似的方法得到（用 $y''$ 代替式中的 $y$ 即得下式），
\begin{eqnarray}
y''(x_{n+1}) + y''(x_{n-1}) - 2 y''(x_n) = h^2 y'''' (x_n) + O(h^4), 
\end{eqnarray}
将（\ref{Numerov1}）代入上式得到，
\begin{eqnarray}
&& - g( x_{n+1} ) y( x_{n+1} ) + s( x_{n+1} ) - g( x_{n-1} ) y( x_{n-1} ) + s( x_{n-1} ) 
\nonumber\\
&& - 2 ( -g(x_n) y(x_n) + s(x_n) ) = h^2 y''''(x_n) + O(h^4), 
\end{eqnarray}
\end{frame}

\begin{frame}{Numerov方法}
将上式代入（\ref{Numerov3}），并整理，即得
\begin{eqnarray}
&& y(x_{n+1}) [ 1 + \frac{h^2}{12} g(x_{n+1}) ] \nonumber\\
&=& 2 y(x_n) [ 1 - \frac{ 5 h^2}{12} g(x_n) ] - y(x_{n-1}) [1 + \frac{h^2}{12} g(x_{n-1}) ]
\nonumber\\
&& + \frac{ h^2 }{12} [ s(x_{n+1}) + 10 s(x_n) + s(x_{n-1}) ] + O(h^6).
\end{eqnarray}
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十章：常微分方程 }
            \begin{itemize}
            \item {\color{blue} 初值问题}
            \item {\color{blue} 边值问题}
            \end{itemize}
		\item 第十一章：偏微分方程
	\end{itemize}
\end{frame}

\begin{frame}{边值问题}
最简单的是两点边值问题，在一个边值处的条件不足以限制唯一解，需要结合另一边的条件才行。
为了直观形象，我们可以再次访问一维无限深势阱薛定谔方程。
\begin{equation}
V(x) = \left\{
\begin{aligned}
&0, ~~~ 0<x<a\\
&\infty, ~~else
\end{aligned}
\right.
\end{equation}
在$[0,a]$区间内，薛定谔方程为
\begin{equation}
- \frac{ \hbar^2 }{2m} \frac{ d^2 }{dx^2} \psi = E \psi，
\end{equation}
相应的本征值和解为
\begin{eqnarray}
E_n &=& \frac{ n^2 \pi^2 \hbar^2 }{ 2 m a^2}, \\
\psi_n &=& \sqrt{\frac{2}{a}} \sin (n \pi x /a  ).
\end{eqnarray}
\end{frame}

\begin{frame}{边值问题}
在数值中，为了简便，我们取$\frac{\hbar^2}{2m} = a = 1$，则薛定谔方程变为
\begin{equation}
- \frac{d^2}{dx^2} \psi = E \psi, 
\end{equation}
其解析解为
\begin{eqnarray}
E_n &=& n^2 \pi^2 , \\
\psi_n &=& \sqrt{2} \sin (n \pi x  ).
\end{eqnarray}

在矩阵本征值一章，我们已经用有限差分将微分方程转化为矩阵问题，求解了这个问题。
在这里，我们将它作为边值问题的典型，用微分方程的方法求解。
在求解时，本征值是未知的，所以不能像初值问题中一样，从左推到右。
\end{frame}

\begin{frame}{打靶法}
我们可以换一个角度看问题，如果已知$\psi(0)=0, \psi'(0)=1$\footnote{这里$\psi'(0)$可以设为不同的值，仅与波函数归一化有关，不影响解的唯一性}，如果已知本征值$E$，则可以用龙格库塔方法从左到右推得$\psi_1$的值。
所以，换一句话说，这样做的话，$\psi_1$是本征值$E$的函数，我们要求的，是使得$\psi_1$为等于$0$的$E$值
\begin{equation}
\psi_1 (E) = 0,
\end{equation}
问题转化为了方程求根问题，每个根对应一个本征值，有了具体的本征值，就可以用龙格库塔方法从左往右推出本征波函数。
如果用牛顿法求解上面的方程，就对应人们比较常说的打靶法。
如果对效率要求不高，也可以用更稳定的二分法。
要做龙格库塔，先把薛定谔方程降阶处理
\begin{eqnarray}
\frac{d}{dx} y^{(0)} &=& y^{(1)}, \\
\frac{d}{dx} y^{(1)} &=& - E y^{(0)}
\end{eqnarray}
然后在程序中定义$\psi_1 (E)$，然后做二分法即可。
\end{frame}

\begin{frame}[fragile]{打靶法}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>
#include"library.h"

double E;//全局变量E

void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];
        dydx[1] = - E* y[0];
}

double psib(double x){

        E=x;//将x的值赋给全局变量E
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{打靶法}
\lstset{firstnumber=last}
\begin{lstlisting}
        int N=10000,n=2;
        double a=0, ya=0, deriv_a=1, b=1;
        int i,j;
        double h = (b-a)/N;
        double *y = new double [n];
        double *ynext = new double [n];
        y[0]=ya; y[1] = deriv_a;
        for(i=0;i<N;i++){
                rk4(a+i*h, h, n, y, ynext, derivative);
                for(j=0;j<n;j++) y[j] = ynext[j];
        }
        return y[0];//返回 psi(b)
}
int main(){

        double eig = bisection(psib, 50, 100, 1E-9);
        cout<<"E/pi/pi="<<eig/M_PI/M_PI<<endl;

        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{应用：三维球对称势场的薛定谔方程}
在球对称势场 $V(r)$ 中，定态薛定谔方程为
\begin{eqnarray}
( - \frac{ \hbar^2 }{ 2m} \nabla^2 + V(r) ) \psi ( \vec{r} ) = E \psi( \vec{r} ),
\end{eqnarray}
假设在球坐标下，体系本征态可以写作
\begin{eqnarray}
\psi( \vec{r} ) = R(r) Y^m_l ( \theta, \phi ),
\end{eqnarray}
则有
\begin{eqnarray}
- \frac{ \hbar^2 }{ 2m } ( \frac{1}{r^2} \frac{ d }{ dr } r^2 \frac{ d }{dr } - \frac{ l(l+1) }{ r^2} ) R(r) + V(r) R(r) = E R(r).
\end{eqnarray}
为了简化上式的形式，可以定义 $R(r) = u(r) / r $，则有
\begin{eqnarray}
- \frac{ \hbar^2 }{ 2m } \frac{ d^2 }{ dr^2 } u(r) + ( V(r) + \frac{ l(l+1) }{ r^2} \frac{ \hbar^2}{2m} ) u(r) = E u(r).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{应用：三维球对称势场的薛定谔方程}
上式还可以进一步化简，定义
\begin{eqnarray}
\kappa &\equiv& \frac{ \sqrt{-2mE} }{ \hbar },  \\
\rho &=& \kappa r, 
\end{eqnarray}
则（43）式变为
\begin{eqnarray}
\frac{ d^2 u }{ d \rho^2 } = [ 1 - \frac{ V(r) }{ E } + \frac{ l(l+1)}{ \rho^2 } ] u(\rho).
\end{eqnarray}
有时这个形式更方便：
\begin{eqnarray}
\frac{ d^2 u }{ d \rho^2 } = [ 1 + \frac{ 2mV( \rho / \kappa) }{ \kappa^2 \hbar^2 } + \frac{ l(l+1)}{ \rho^2 } ] u(\rho).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{应用：三维球对称势场的薛定谔方程}
解这个问题的算法如下：
\begin{shaded*}
\begin{itemize}
\item [1] 确定边界条件。
在无穷远处，有 $u(\infty) = 0$。
由于 $R(r)$ 在 $r=0$ 处为有限值，所以 $u(r) = r R(r) $ 在 $r=0$ 处有 $u(0) = 0$，而 $u'(0)$ 可以任取，只会影响 $u(r)$ 的归一性。
\item [2] 选取一个足够大的值 $r_{max}$，将 $[0, r_{max}]$ 等分为 $n$ 份。
设定初始值 $u(r_{max}) = 1E-5, u'(r_{max}) = 0$，用 numerov 方法，从右往左推知各节点的函数值。
若未遇到极值点，则在 $r_{max}/10$ 处停止；若遇到极值点，则在极值点处停止。
设停止处为 $r_m$，从右往左推得该处函数值为 $u_{<} ( r_m )$。
\item [3] 用 numerov 方法，从 $u(0)=0, u'(0) = 1$ 推得 $r_m$ 处函数值 $u_{>}(r_m)$。
\item [4] 将 $r_m$ 右侧所有函数值乘以 $ u_< (r_m) / u_> (r_m) $，即使得函数连续。
\item [5] 调节本征值 $\rho_0$，使得 $ u (r_m) - u( r_m - h ) = u( r_m + h) - u( r_m )$，即使得一阶导数连续。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{氢原子能级}
氢原子的势场为
\begin{eqnarray}
V(r) = - \frac{ e^2 }{ 4 \pi \epsilon_0 r },
\end{eqnarray}
其中 $e = 1.6 \times 10^{-19} C$。
将 $V(r)$ 代入（44），得到
\begin{eqnarray}
\frac{ d^2 u }{ d \rho^2 } = [ 1 - \frac{ \rho_0 }{ \rho } + \frac{ l(l+1)}{ \rho^2 } ] u,
\end{eqnarray}
其中
\begin{eqnarray}
\rho_0 &=& \frac{ me^2 }{ 2 \pi \epsilon_0 \hbar^2 \kappa}.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{氢原子能级}
能量与 $\rho_0$ 的换算关系为
\begin{eqnarray}
 E = - [ \frac{ m }{ 2 \hbar^2 } ( \frac{ e^2 }{ 4 \pi \epsilon_0 } )^2 ] \frac{ 4 }{ \rho^2_0 }.
\end{eqnarray}

而理论上氢原子的本征能量为
\begin{eqnarray}
 E = - [ \frac{ m }{ 2 \hbar^2 } ( \frac{ e^2 }{ 4 \pi \epsilon_0 } )^2 ] \frac{ 1 }{ n^2 },
\end{eqnarray}
其中 $n$ 为主量子数，轨道角动量量子数 $l = n-1, n-2, ..., 0$。

所以，$\rho_0$ 的理论值为 $\rho^{th}_0 = 2, 4, 6, 8, ...$。
表中展示了计算出来的 $\rho^{cal}_0$ 的值，与理论值相符。
\end{frame}

\begin{frame}[fragile]{氢原子能级}
\begin{table}
\caption{ 氢原子能谱对应的 $\rho_0$ 计算值 $\rho^{cal}_0$，$\rho_0$ 的定义请见文中。
在计算中，我们选取 $\rho_{max} = 50$，将 [0, 50] 分为 10000 个区间。
}
\centering
\begin{tabular}{ccccc}
\hline \hline
	$l=0$		&	$l=1$		&	$l=2$		& $l=3$		&	$l=4$		\\
\hline
	2.00018	&				&				&				&				\\
	4.00063	&	3.99989	&				&				&				\\
	6.0025	&	6.0004	&	5.99992	&				&				\\
	8.00451	&	8.00154	&	8.00154	&	8.00054	&				\\
	10.006	&	10.0014	&	10.0014	&	10.0014	&	10.0013	\\
	12.0077	&	12.0013	&	12.0013	&	12.0013	&	12.0013	\\
	14.0098	&	14.0012	&	14.0013	&	14.0013	&	14.0013	\\
	16.0121	&	16.0011	&	16.0012	&	16.0012	&	16.0012	\\
	18.0147	&	18.001		&	18.0012	&	18.0012	&	18.0012	\\
	20.0175	&	20.0009	&	20.0011	&	20.0011	&	20.0011	\\
\hline \hline
\end{tabular}
\end{table}
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十章：常微分方程 }
            \begin{itemize}
            \item {\color{blue} 初值问题}
            \item {\color{blue} 边值问题}
            \end{itemize}
		\item {\color{blue}第十一章：偏微分方程 }
            \begin{itemize}
            \item 扩散方程
            \item 泊松方程
            \item 一维弦振动方程
            \end{itemize}
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{扩散方程}
在数学物理方法的课程中，我们学过热量扩散方程,
\begin{equation}
\nabla^2 T  = D \frac{ \partial T }{ \partial t },
\end{equation}
如果稍微重新定义$T$，即可把微分方程写成如下形式,
\begin{equation}
\nabla^2 u = \frac{ \partial u }{ \partial t }，
\end{equation}
一维情况下，扩散方程为，
\begin{equation}
\frac{ \partial^2 u }{ \partial x^2 } = D \frac{ \partial u }{ \partial t },
\end{equation}
$u(x,t)$在$x$轴上的取值范围是$[x_a,x_b]$，在$t$轴上的取值范围是$[t_a,t_b]$，在初始时刻,
\begin{equation}
u(x,t=t_a) = g(x),
\end{equation}
在边界上，
\begin{eqnarray}
u(x_a,t) &=& a(t), \\
u(x_b,t) &=& b(t).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{显式递推公式}
利用有限差分，可以近似处理偏导数，
\begin{eqnarray}
\frac{ \partial^2 u }{ \partial x^2 } &=& \frac{1}{(\Delta x)^2}( u(x+ \Delta x,t) + u(x- \Delta x,t) - 2u(x,t) ) + O((\Delta x)^2),
\\
\frac{ \partial u }{ \partial t } &=& \frac{1}{\Delta t}( u(x,t+ \Delta t) - u(x,t) ) + O(\Delta t).
\end{eqnarray}
带入一维扩散方程，得到
\begin{equation}
\frac{1}{(\Delta x)^2}( u(x+ \Delta x,t) + u(x- \Delta x,t) - 2u(x,t) )
\approx \frac{1}{\Delta t}( u(x,t+ \Delta t) - u(x,t) ),
\end{equation}
误差为$O((\Delta x)^2) + O(\Delta t)$。

记$\alpha = \frac{ \Delta t }{ (\Delta x)^2 }$，则方程变为
\begin{equation}
u(x,t+ \Delta t) \approx \alpha u(x+ \Delta x,t) + (1- 2\alpha)u(x,t)+ \alpha u(x- \Delta x,t),
\end{equation}
\end{frame}

\begin{frame}[fragile]{显式递推公式}
这意味着，如果我们知道了$t$时刻的情况，就可以往前推断$\Delta t$之后的情况。
若$u(x,t)$在$x$轴上的取值范围是$[x_a,x_b]$，在$t$轴上的取值范围是$[t_a,t_b]$，我们把这个定义域用四方网格划分，格点在$x$方向上间距为$\Delta x$，在$t$方向上的间距为$\Delta t$，记
\begin{eqnarray}
x_i &=& x_a + i*\Delta x, \\
t_j &=& t_a + j*\Delta t,
\end{eqnarray}
格点上的$u$值标记为
\begin{equation}
u_{i,j} = u(x_i, t_j) \approx u(x_a + i*\Delta x, t_a + j * \Delta t ),
\end{equation}
那么，数值递推式为
\begin{equation}
u_{i,j+1} = \alpha u_{i-1,j} + (1- 2\alpha)u_{i,j} + \alpha u_{i+1,j},
\end{equation}
根据前人的推断[cite Golub]，这个公式在如下条件下是稳定的，
\begin{equation}
\Delta t / (\Delta x)^2 \leq 1/2,
\end{equation}
所以如果$\Delta t$太大，结果可能是错误的。
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
在上面的推导中，我们计算$u$对时间的偏导数时，用了公式
\begin{equation}
\frac{ \partial u }{ \partial t } = \frac{1}{\Delta t}( u(x,t+ \Delta t) - u(x,t) ) + O(\Delta t).
\end{equation}
我们也可以用这个公式
\begin{equation}
\frac{ \partial u }{ \partial t } = \frac{1}{\Delta t}( u(x,t) - u(x,t- \Delta t) ) + O(\Delta t).
\end{equation}
那么，经过类似的推导，并整理后，得到数值递推式为
\begin{equation}
u_{i,j-1} = -\alpha u_{i-1,j} + (1+ 2\alpha)u_{i,j} - \alpha u_{i+1,j},
\end{equation}
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
如果$u_{*,j-1}$都已经已知，我们可以得到关于$u_{*,j}$的$N_x$个线性方程，$N_x$为$x$轴上$(x_a,x_b)$之间的格点数。
\begin{equation}
\begin{bmatrix}
1+2\alpha & -\alpha \\
-\alpha   & 1+2\alpha & -\alpha \\
          & -\alpha   & 1+2\alpha & -\alpha \\
          &			  &	\ddots	  & \ddots  & \ddots \\
          &			  &			  & -\alpha & 1+2\alpha
\end{bmatrix}
\begin{bmatrix}
u_{1,j} \\
u_{2,j} \\
\vdots \\
u_{N_x,j}
\end{bmatrix}
=
\begin{bmatrix}
u_{1,j-1} + \alpha a(t_j)\\
u_{2,j-1} \\
\vdots\\
u_{N_x,j-1} + \alpha b(t_j)
\end{bmatrix}
\end{equation}
于是我们用三对角线性方程组的数值解法，求解上面的方程组，就得到了$u_{*,j}$，因为上式没有写出$u_{*,j}$的显式表达式，所以这种方法叫做隐式递推。
前文中提到过，解三对角线性方程组只需要$O(N_x)$次操作，所以是比较有效的。
根据前人的证明，这个算法是稳定的。
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
下面是隐式递推的代码，征用了前面章节编写的三对角方程组求解函数TridiagLinear。
\begin{lstlisting}
/*
 * diffusion1d(...)             solves a 1-dimensional diffusion equation
 *              \partial^2 u / \partial x^2 = \partial u / \partial t
 *              and prints u(*,t_j) into "diffusion1d.txt"
 * double xa
 * double xb                    x zone [xa, xb]
 * double (*g_initial)(...)     u(x, t=0)
 * double ta
 * double tb                    t zone [ta, tb]
 * double (*boundary_a)(...)    u(xa, t)
 * double (*boundary_b)(...)    u(xb, t)
 * int Nx                       grid point number in (xa,xb)
 * int Nt                       grid point number in (ta,tb)
 */
void diffusion1d( double xa, double xb, double (*g_initial)(double x), double ta, double tb, double (*boundary_a)(double t), double (*boundary_b)(double t), int Nx, int Nt){
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
\lstset{firstnumber = last}
\begin{lstlisting}
        int i, j;
        double dx, dt, alpha;
        double *u = new double [Nx];
        double *unext = new double [Nt];
        double *a = new double [Nx];// tridiagonal matrix
        double *b = new double [Nx];
        double *c = new double [Nx];

        dx = (xb-xa)/(Nx+1);// [xa, xb] is divided into Nx+1 small zones
        dt = (tb-ta)/Nt;// [ta, tb] is divided into Nt steps
        alpha = dt/dx/dx;
        for(i=0;i<Nx;i++){// initialize tridiagonal matrix
                a[i] = -alpha;
                b[i] = 1 + 2*alpha;
                c[i] = -alpha;
        }
        for(i=0;i<Nx;i++){//initial u(x, t=0)
                u[i] = g_initial(xa + (i+1)*dx);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
\lstset{firstnumber = last}
\begin{lstlisting}
        ofstream fp("diffusion1d.txt");
        for(j=1;j<=Nt;j++){
                fp<< boundary_a( ta+ (j-1)*dt)<<"\t";//output u(*,j-1) into file
                for(i=0;i<Nx;i++)fp<<u[i]<<"\t";
                fp<< boundary_b( ta+ (j-1)*dt)<<"\t";
                fp<<endl;
                u[0] += alpha * boundary_a( ta + j*dt);
                u[Nx-1] += alpha * boundary_b( ta + j*dt);
                TridiagLinear( Nx, a, b, c, unext, u);
                for(i=0;i<Nx;i++) u[i] = unext[i];
        }
        fp.close();
        delete [] c;
        delete [] b;
        delete [] a;
        delete [] u;
        delete [] unext;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
自定义初始条件，边界条件以后，就可以调用上面的代码计算1维扩散问题。
下面的代码做一个热量传导问题，
\begin{lstlisting}
头文件。。。
函数diffusion1d(...)

double Tleft(double t){// temperature of left boundary
        return 100;
}

double Tright(double t){// temperature of right boundary
        return 0;
}

double Tinitial(double x){
        return 0;
}

int main(){
        diffusion1d(0, 1, Tinitial, 0, 0.5, Tleft, Tright, 1000, 1000);
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{隐式递推公式}
一根铁棒上各处温度都是0摄氏度，右端浸泡在0摄氏度的水里，左端与100摄氏度的沸水接触，随着时间的过去，热量会源源不断地从铁棒左端流向右端，铁棒上的温度分布会趋于均匀。

运行上面的代码以后，生成diffusion1d.txt文件，然后用gnuplot，可以画出热量图，如图\ref{fig:diffusion1d}
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{diffusion1d}
\caption{一维热扩散问题：铁棒上的温度演化图。铁棒左端恒为100摄氏度，右端恒为0摄氏度，初始时铁棒上温度为0。}
\label{fig:diffusion1d}
\end{figure}
\end{frame}

\begin{frame}[fragile]{二维泊松方程}
二维泊松方程为
\begin{equation}
\nabla^2 u = \rho,
\end{equation}
静电势的方程就是这个形式。
这里我们求解二维情形，三维求解的道理也一样，只是计算量比较大，算一次要等很久。

二维情况下，泊松方程变为
\begin{equation}
\frac{\partial^2 }{\partial x} u + \frac{\partial^2}{\partial y} u = \rho.
\end{equation}
用有限差分公式代替二阶偏导数，得到
\begin{equation}
\frac{ u_{i+1,j} + u_{i-1,j} - 2 u_{i,j} }{(\Delta x)^2} + O((\Delta x)^2)
+ \frac{ u_{i,j+1} + u_{i,j-1} - 2 u_{i,j} }{(\Delta y)^2} + O((\Delta y)^2) = \rho_{i,j}
\end{equation}
其中，$u_{i,j} = u(x_a + i\Delta x, y_a + j\Delta y)$.
\end{frame}

\begin{frame}[fragile]{二维泊松方程}
那么，如果取$\Delta x = \Delta y = h$，在2阶精度下，数值递推式为
\begin{equation}
u_{i-1,j} + u_{i+1,j} + u_{i,j+1}, u_{i,j-1} - 4 u_{i,j} = h^2 \rho_{i,j}.
\end{equation}
写成矩阵会得到一个关于$u_{i,j}$的线性方程组，系数矩阵是一个稀疏矩阵，求得线性方程组的解，就得到任意格点的$u_{i,j}$，就认为得到了问题的数值解。

如果用高斯消元法，或者LU分解，会得到大量的非零矩阵元，使整个问题变得更加复杂。
所以这种问题，更适合用迭代法。
如果用Jacobi方法，即
\begin{eqnarray}
A \vec{x} &=& \vec{b}, \\
A &=& U + D + L, \\
\vec{x}_{j+1} &=& D^{-1}( \vec{b} - (U+L) \vec{x}_j ).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{二维泊松方程}
给每个格点上的$u_{i,j}$任意赋一个初始值，然后用Jacobi方法进行迭代，实际上就相当于
\begin{equation}
u^{(k+1)}_{i,j} = \frac{1}{4}(u^{(k)}_{i-1,j} + u^{(k)}_{i+1,j} + u^{(k)}_{i,j+1}, u^{(k)}_{i,j-1} - h^2 \rho_{i,j}).
\end{equation}
如果在计算等号左边的$u^{(k+1)}_{i,j}$时，右边的式子中有些项已经做了第$k+1$轮迭代，比如，已经有$u^{(k+1)}_{i-1,j}$，那么可以不用$u^{(k)}_{i-1,j}$，而用$u^{(k+1)}_{i-1,j}$，这样做就是高斯-席德尔方法。

我们可以用高斯席德尔方法，写一个通用的二维四方网格扩散方程求解函数，然后做一个简单的例子：匀强电场。
\end{frame}

\begin{frame}[fragile]{二维泊松方程}
\begin{lstlisting}
void poisson2d(int n, double **u, double xa, double ya, double h, double precision,
                double (* boundary_x1)(double x), double (* boundary_x2)(double x),
                double (* boundary_y1)(double y), double (* boundary_y2)(double y),
                double (* rho)(double x, double y) ){

        int i,j;
        double temp, change = precision +1;

        for(i=0;i<n;i++) u[i][0] = boundary_x1( xa + i*h );
        for(i=0;i<n;i++) u[i][n-1] = boundary_x2( xa + i*h );
        for(j=0;j<n;j++) u[0][j] = boundary_y1( ya + i*h );
        for(j=0;j<n;j++) u[n-1][j] = boundary_y2( ya + i*h );

        for(i=1;i<n-1;i++)
                for(j=1;j<n-1;j++)
                        u[i][j]=0;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{二维泊松方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        while( change > precision ){

                change = 0;
                for(i=1;i<n-1;i++){
                        for(j=1;j<n-1;j++){
                                temp = u[i][j];
                                u[i][j] = 0.25 * (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] - h * h * rho(xa+i*h, ya+i*h) );
                                change += fabs( u[i][j] - temp );
                        }
                }
        }
        ofstream fp("poisson2d.txt");
        for(j=0;j<n;j++){
                for(i=0;i<n;i++){
                        fp<<u[i][j]<<"\t";
                }
                fp<<endl;
        }
        fp.close();
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{匀强电场}
要做匀强电场，需要设置恰当的边界条件，
\begin{lstlisting}
double x1(double x){
        return x;
}
double x2(double x){
        return x;
}
double y1(double y){
        return 0;
}
double y2(double y){
        return 1;
}
double rho(double x, double y){
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{匀强电场}
\lstset{firstnumber = last}
\begin{lstlisting}
int main(){

        int i,n=100;
        double xa=0, ya=0, h=1.0/n;

        double **u = new double *[n];
        for(i=0;i<n;i++) u[i] = new double [n];

        poisson2d( n, u, xa, ya, h, 0.001, x1, x2, y1, y2, rho);
        return 0;
}
\end{lstlisting}
这里设置左边界为电势为$0$的极板，右边界为电势为$1$的极板，上下边界电势均匀变化，所以会得到匀强电场。
\end{frame}

\begin{frame}[fragile]{匀强电场}
得到数据以后，用gnuplot作图
\begin{lstlisting}
set term png
set output "poisson2d.png"
unset border
set xtics ("0.0" 0, "0.2" 20, "0.4" 40, "0.6" 60, "0.8" 80, "1.0" 99)
set ytics ("0.0" 0, "0.2" 20, "0.4" 40, "0.6" 60, "0.8" 80, "1.0" 99)
plot "poisson2d.txt" matrix with image
\end{lstlisting}
就会得到二维电势图\ref{fig:poisson2d}。
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{poisson2d}
%\caption{二维泊松方程：匀强电场。左侧极板电势为0，右侧极板电势为1。}
\label{fig:poisson2d}
\end{figure}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
一维波动方程（比如绳波）的方程为
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\label{eqn:wave1d}
\end{eqnarray}
如果初始时1m长的绳子静止，绳子一端做简谐振动，另一端固定在墙上，那么初始条件和边界条件分别为
\begin{eqnarray}
u(x,t=0) &=& 0, \\
u(0,t) &=& sin(t), \\
u(1,t) &=& 0.
\end{eqnarray}
还需要一个初始条件，描述初始时刻绳子上各个质点的运动速度，
\begin{equation}
\frac{\partial}{\partial t} u |_{t=0} = 0.
\end{equation}
如果把微分方程用有限差分代替，即
\begin{eqnarray}
\frac{ u_{i+1,j} + u_{i-1,j} - 2 u_{i,j} }{(\Delta x)^2} + O((\Delta x)^2)
= \frac{ u_{i,j+1} + u_{i,j-1} - 2 u_{i,j} }{(\Delta t)^2} + O((\Delta t)^2),
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
如果取$\alpha = \frac{(\Delta t)^2}{(\Delta x)^2}$，显式的递推式为[cite Morten]
\begin{equation}
u_{i,j+1} = (2-2\alpha) u_{i,j} + \alpha u_{i+1,j} + \alpha u_{i-1,j} - u_{i,j-1}.
\end{equation}
只要知道$u_{*,j-1}$与$u_{*,j}$，就可以递推得出$u_{*,j+1}$。

第一步前推时，我们只知道$u_{*,0}$，不知道$u_{*,-1}$，那么可以利用初始条件
\begin{equation}
\frac{\partial}{\partial t} u |_{t=0} = 0.
\end{equation}
翻译成差分，可以得到
\begin{equation}
u_{i,-1} = u_{i,1},
\end{equation}
带入递推式，得到
\begin{equation}
u_{i,1} = (1-\alpha) u_{i,0} + \alpha \frac{u_{i+1,0} + u_{i-1,0}}{2}.
\end{equation}
这样就可以完成前推了。
\end{frame}

\begin{frame}[fragile]{1维波动方程}
我还不清楚这个显式格式的稳定性，所以谨慎一点，将$\Delta t$取得小一点看看。
下面的代码中是一个函数，在给定初始条件，边界条件下，计算一维波动方程(\ref{eqn:wave1d})，并将各个时间节点的所有$u_{*,j}$作为一行输出到给定文件"file"中。
\begin{lstlisting}
void wave1d(    double xa, int nx, double dx,
                double (*boundary1)(double t), double (*boundary2)(double t),
                int nt, double dt, double (* u_initial)(double x),
                string file ){

        int i,j;
        double alpha = dt*dt/dx/dx;

        double *u_old = new double [nx]; // u_{i-1}
        double *u = new double [nx]; // u_i
        double *u_new = new double [nx]; // u_{i+1}

        u_old[0] = boundary1(0);
        u_old[nx-1] = boundary2(0);
        for(i=1;i<nx-1;i++) u_old[i] = u_initial( xa + i*dx);
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
\lstset{firstnumber=last}
\begin{lstlisting}
        u[0] = boundary1(dt);
        u[nx-1] = boundary2(dt);
        for(i=1;i<nx-1;i++) u[i] = (1-alpha)*u_old[i] + alpha/2*( u_old[i+1] + u_old[i-1] );

        ofstream fp(file);
        if(!fp){
                cout<<"error: failed to open "<<file<<endl;
                exit(1);
        }
        for(i=0;i<nx;i++)fp<<u_old[i]<<"\t"; fp<<endl;
        for(i=0;i<nx;i++)fp<<u[i]<<"\t"; fp<<endl;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
\lstset{firstnumber=last}
\begin{lstlisting}
        for(j=2;j<=nt;j++){
                u_new[0] = boundary1( j*dt);
                u_new[nx-1] = boundary2( j*dt);
                for(i=1;i<nx-1;i++)
                        u_new[i] = (2-2*alpha)*u[i] + alpha * u[i+1] + alpha* u[i-1] - u_old[i];
                for(i=0;i<nx;i++)fp<<u_new[i]<<"\t"; fp<<endl;

                for(i=0;i<nx;i++){
                        u_old[i] = u[i];
                        u[i] = u_new[i];
                }
        }
        fp.close();
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
我们可以设置初始时绳子静止，左端点开始做简谐振动，右端点固定。
那么简谐波会从左往右传播，在右端点处反弹，反弹的波从右往左传播，与从左往右的波相干涉，形成驻波。
生成数据文件的代码如下：
\begin{lstlisting}

double u_initial(double x){
        return 0;
}

double boundary1(double t){
        return sin(10*t);
}

double boundary2(double t){
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
\lstset{firstnumber=last}
\begin{lstlisting}
int main(){

        int i,j,nx=300;
        double xa=0, xb=3;
        double dx = (xb-xa)/nx;
        int nt=500;
        double ta=0, tb=5;
        double dt = (tb-ta)/nt;
        double **matrix;

        string file = "wave1d.txt";

        wave1d( 0, nx, dx, boundary1, boundary2,
                nt, dt, u_initial, file);
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
\lstset{firstnumber=last}
\begin{lstlisting}
        // transpose the data matrix
        matrix = new double *[nt+1];
        for(i=0;i<nt+1;i++) matrix[i] = new double [nx];

        ifstream fp1( file );
        for(i=0;i<nt+1;i++)
                for(j=0;j<nx;j++)
                        fp1>> matrix[i][j];
        fp1.close();

        ofstream fp2( file );
        for(j=0;j<nx;j++){
                fp2<< xa + j*dx <<"\t";
                for(i=0;i<nt+1;i++)
                        fp2<< matrix[i][j]<<"\t";
                fp2<<endl;
        }
        fp2.close();

        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
代码中最后一段是将数据文件的行和列做转置，方便用gnuplot画成不同的桢，然后生成动图gif文件。
gnuplot代码如下
\begin{lstlisting}
set term gif animate delay 1
set output "wave1d.gif"
set xrange [0:3]
set yrange [-2:2.5]
set xlabel "x(m)"
set ylabel "y(m)"
do for[t=2:502:1]{
        plot "wave1d.txt" u 1:t w l lw 3 t sprintf("string wave: v=1m/s, t=%f s", (t-2)*0.01)
}
set output
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
2维波动方程如下
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u + \frac{ \partial^2 }{ \partial y^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\label{eqn:wave2d}
\end{eqnarray}
初始条件和边界条件为
\begin{eqnarray}
u(x,y,t=0) &=& g(x,y), \\
\frac{\partial}{\partial t} u(x,y,t) |_{t=0} &=& 0, \\
u(x=x_a, y, t) &=& Bxa(y,t), \\
u(x=x_b, y, t) &=& Bxb(y,t), \\
u(x, y=y_a, t) &=& Bya(x,t), \\
u(x, y=y_b, t) &=& Byb(x,t).
\end{eqnarray}
如果取$\Delta x = \Delta y = h$，则用有限差分代替2阶偏微分，可以得到
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\begin{eqnarray}
&&\frac{ u^k_{i+1,j} + u^k_{i-1,j} - 2 u^k_{i,j} }{h^2} + O(h^2)
+ \frac{ u^k_{i,j+1} + u^k_{i,j-1} - 2 u^k_{i,j} }{h^2} + O(h^2) \nonumber\\
&=& \frac{ u^{k+1}_{i,j} + u^{k-1}_{i,j} - 2 u^k_{i,j} }{(\Delta t)^2} + O((\Delta t)^2),
\end{eqnarray}
其中$u^k_{i,j} = u(x_a + i*h, y_a + j*h, t_a + k*\Delta t)$。
记$\alpha = (\Delta t)^2 / h^2$，得到显式递推式
\begin{equation}
u^{k+1}_{i,j} = (2-4\alpha)u^k_{i,j} - u^{k-1}_{i,j} + \alpha ( u^k_{i-1,j} + u^k_{i+1,j} + u^k_{i,j-1} + u^k_{i,j+1} ).
\end{equation}
在初始时，利用初始条件$\frac{ \partial }{ \partial t} u |_{t=0} = 0$，可以得到 $u^{-1}_{i,j} = u^1_{i,j}$，得到递推公式的第一步
\begin{equation}
u^1_{i,j} = (1-2\alpha)u^0_{i,j} + \frac{\alpha}{2} ( u^0_{i-1,j} + u^0_{i+1,j} + u^0_{i,j-1} + u^0_{i,j+1} ).
\end{equation}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
所以很容易将前面的1维波动代码改写为2维波动代码。
下面的代码实现二维平面波干涉，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include<string>

double source(double t){
        return sin(2*M_PI*t);
}

void wave2d(    double xa, double ya, int n, double h,
                double (*Bxa)( double y, double t ), double (*Bxb)( double y, double t ),
                double (*Bya)( double x, double t ), double (*Byb)( double x, double t ),
                int nt, double dt, double (* u_initial)(double x, double y),
                string file ){
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        int i,j,k;
        double alpha = dt*dt/h/h;
        int ct = n/2;

        ofstream gnu("wave2d_interference.gnu");
        gnu<<"set term gif animate delay 5\n";
        gnu<<"set output \"wave2d_interference.gif\"\n";
        gnu<<"set xrange [0:10]\n";
        gnu<<"set yrange [0:10]\n";
        gnu<<"set zrange [-0.5:1]\n";
        gnu<<"unset border\n";
        gnu<<"unset xtics\n";
        gnu<<"unset ytics\n";
        gnu<<"unset ztics\n";
        gnu<<"set view 45, 30, 1, 1\n";
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        double **u_old = new double *[n];
        for(i=0;i<n;i++) u_old[i] = new double [n];
        double **u = new double *[n];
        for(i=0;i<n;i++) u[i] = new double [n];
        double **u_new = new double *[n];
        for(i=0;i<n;i++) u_new[i] = new double [n];

        // calculate u^0_{i,j}
        for(i=0;i<n;i++){
                u_old[0][i] = Bxa( ya + i*h , 0 );
                u_old[n-1][i] = Bxb( ya + i*h, 0 );
                u_old[i][0] = Bya( xa + i*h, 0 );
                u_old[i][n-1] = Byb( xa + i*h, 0 );
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        u_old[ct-15][ct] = source(0);
        u_old[ct+15][ct] = source(0);
        for(i=1;i<n-1;i++){
                for(j=1;j<n-1;j++){
                        if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
                        u_old[i][j] = u_initial( xa + i*h, ya + i*h );
                }
        }
        // output u^0_{i,j}
        string file_frame = file + "000.txt";
        ofstream fp(file_frame);
        if(!fp){
                cout<<"error: failed to open "<<file_frame<<endl;
                exit(1);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u_old[i][j]<<endl;
                }
                fp<<endl;
        }
        fp.close();

        gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", 0)\n";


        // calculate u^1_{i,j}
        for(i=0;i<n;i++){
                u[0][i] = Bxa( ya + i*h , dt );
                u[n-1][i] = Bxb( ya + i*h, dt );
                u[i][0] = Bya( xa + i*h, dt );
                u[i][n-1] = Byb( xa + i*h, dt );
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        u[ct-15][ct] = source(dt);
        u[ct+15][ct] = source(dt);
        for(i=1;i<n-1;i++){
                for(j=1;j<n-1;j++){
                        if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
                        u[i][j] = (1-2*alpha)*u_old[i][j] + alpha/2*( u_old[i-1][j] + u_old[i+1][j] + u_old[i][j+1] + u_old[i][j-1] );
                }
        }

        // output u^1_{i,j}
        file_frame = file + "001.txt";
        fp.open(file_frame);
        if(!fp){
                cout<<"error: failed to open "<<file_frame<<endl;
                exit(1);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u[i][j]<<endl;
                }
                fp<<endl;
        }
        fp.close();

        gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", "<<dt<<")\n";
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(k=2;k<=nt;k++){
                // calculate u^{k}_{i,j}
                for(i=0;i<n;i++){
                        u_new[0][i] = Bxa( ya + i*h , k*dt );
                        u_new[n-1][i] = Bxb( ya + i*h, k*dt );
                        u_new[i][0] = Bya( xa + i*h, k*dt );
                        u_new[i][n-1] = Byb( xa + i*h, k*dt );
                }
                u_new[ct-15][ct] = source(k*dt);
                u_new[ct+15][ct] = source(k*dt);
                for(i=1;i<n-1;i++){
                        for(j=1;j<n-1;j++){
                                if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
                                u_new[i][j] = (2-4*alpha)*u[i][j] - u_old[i][j] + alpha*( u[i-1][j] + u[i+1][j] + u[i][j+1] + u[i][j-1] );
                        }
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
                for(i=0;i<n;i++){
                        for(j=0;j<n;j++){
                                u_old[i][j] = u[i][j];
                                u[i][j] = u_new[i][j];
                        }
                }
                // output u^k_{i,j}
                file_frame = file + (char)(48+(k/100)%10) + (char)(48+ (k/10)%10) + (char)(48+ (k%10)) + ".txt";
                fp.open(file_frame);
                if(!fp){
                        cout<<"error: failed to open "<<file_frame<<endl;
                        exit(1);
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
                for(i=0;i<n;i++){
                        for(j=0;j<n;j++){
                                fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u_new[i][j]<<endl;
                        }
                        fp<<endl;
                }
                fp.close();
                gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", "<< k*dt <<")\n";
        }
        gnu<<"set output"<<endl;
        gnu.close();
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
double u_initial(double x, double y){
        return 0;
}

double Bxa(double y,double t){
        return 0;
}
double Bxb(double y, double t){
        return 0;
}
double Bya(double x, double t){
        return 0;
}
double Byb(double x, double t){
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维波动方程}
\lstset{firstnumber = last}
\begin{lstlisting}
int main(){

        int i,j,n=300;
        double xa=0, xb=10;
        double ya=xa, yb=xb;
        double h = (xb-xa)/n;
        int nt=300;
        double ta=0, tb=5;
        double dt = (tb-ta)/nt;
        cout<<"dt="<<dt<<endl;

        string file = "wave2d/interference";
        wave2d( xa, ya, n, h,
                Bxa, Bxb, Bya, Byb,
                nt, dt, u_initial, file );

        return 0;
}
\end{lstlisting}
这个代码既做微分方程的计算，生成数据文件，也生成gnuplot脚本，运行生成的gnuplot脚本，会生成2维平面波干涉的gif动图。
\end{frame}

\begin{frame}{作业}
1. $x(t)$满足微分方程
\begin{equation}
\frac{dx}{dt} = -x,
\end{equation}
且有初始值$x(0)=1$，试用欧拉法求得$[0,5]$范围内的解$x(t)$，将数值解与解析解$e^{-x}$画在同一张图上。

2. 对于平抛运动，其运动方程为
\begin{eqnarray}
\frac{d}{dt} x &=& v, \\
\frac{d^2}{dt^2} y &=& -g,
\end{eqnarray}
不妨设$v=1,g=10$，初始时刻$x=0,y=0$，
\begin{itemize}
\item [1)] 试用龙格库塔方法，求取质点$t=10s$时的坐标。
\item [*2)] 稍微修改程序，求得$y=-h=-5$时的$x$值。
\end{itemize}
\end{frame}

\begin{frame}{作业}
3. 一维方势井可以加一个微扰，使得
\begin{equation}
V(x) = \left\{
\begin{aligned}
&0,~~ &-1<x<0\\
&0.1,~~ &0<x<1\\
&\infty,~~ &else
\end{aligned}
\right.
\end{equation}
取$\hbar = m = 1$, 试求解基态能量。
\\
提示：一阶微扰近似得到$E_0 = \frac{\pi^2}{8} + 0.05$，这个结果不是严格解，但应该有几位有效数字的精度，所以计算结果不可能距离它很远。
\end{frame}

\begin{frame}{作业}
*4. 一维谐振子势为
\begin{equation}
V(x) = \frac{1}{2} k x^2,
\end{equation}
为了叙述简便，不妨取$\hbar = m = k = 1$, 则薛定谔方程为
\begin{eqnarray}
- \frac{d^2}{dx^2} \psi + x^2 \psi &=& 2 E \psi, \\
                                 E &=& (n+\frac{1}{2}).
\end{eqnarray}
可以取$[-5,5]$区间，设置左侧边界条件为$\psi(-5)=0, \psi'(-5)=0.01$，右边界条件为$\psi(5)=0$。求解三个最低本征值与本征波函数。

\end{frame}

\begin{frame}{作业}

5.　试讨论非齐次(有源)热传导问题的数值解法(写出算法思路和核心公式，表述清楚即可)：
\begin{equation}
\left\{
\begin{aligned}
& u_t = a^2 u_{xx} + f(x,t)      & (0 \leq x \leq a, t>0), \\
& u(0,t) = \mu_1 (t); u(a,t) = \mu_2 (t)  & (t \geq 0), \\
& u(x,0) = \varphi(x);           &  ( 0 \leq x \leq a)
\end{aligned}
\right.
\end{equation}

6. 绳子长为3m，两端固定，初始时刻，绳子上各处质点的横向位移为
\begin{equation}
g(x,t=0) = \left\{
\begin{aligned}
& \sin (4 \pi x), & 1<x<1.25, \\
& 0, &else
\end{aligned}
\right.
\end{equation}
绳子的波动方程为
\begin{equation}
\frac{\partial^2}{\partial x^2} u = \frac{\partial^2}{\partial t^2} u. 
\end{equation}
试求5秒之后绳子的形状。
\end{frame}

\end{document}